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Photon-number resolving detectors are a fundamental building-block of optical quantum informa- 
tion processing protocols. A loop detector, combined with appropriate statistical processing, can be 
used to convert a binary on/off photon counter into a photon- number-resolving detector. Here we 
describe the idea of a signature of photon-counts, which may be used to more robustly reconstruct 
the photon number distribution of a quantum state. The methodology is applied experimentally 
in a 9-port loop detector operating at a telecommunications wavelength and compared directly to 
the approach whereby only the number of photon-counts is used to reconstruct the input distribu- 
■ tion. The signature approach is shown to be more robust against calibration errors, exhibit reduced 

f"^ ' statistical uncertainty, and reduced reliance on a-priori assumptions about the input state. 
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I. INTRODUCTION 
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One of the most promising applications of quantum science lies in the field of quantum information processing [l[ . 
Quantum optics provides an extremely convenient test bed for quantum information ideas because quantum effects 
dominate even at room temperatures and non-classical states are relatively easy to achieve in the laboratory. This 
C^h 1 has been most impressively illustrated with recent small-scale experimental demonstrations of quantum algorithms 

ai. 

Photon-number resolving detectors, as well as the ability to characterise quantum optical states using such detectors, 
are a key element in the further development of quantum optical information processing systems Due to the 

low energies carried by photons, contemporary detectors are incapable of resolving photon number with high fidelity. 
Typically linear PIN photodiodes used as direct detectors of the optical field are limited to mesoscopic photon- 
number regimes @|. Recent results show promise at telecommunications wavelengths using liquid-He cooled diodes 
0, bolometric devices [1] and super-conducting nano-structures 0] although such solutions are not yet widespread 

^ i , due to their inherent cost and/or cryogenic cooling complexities. 

Consequently, use is commonly made of single photon detector modules (SPDMs), based on avalanche photo-diodes 
&\ , (APDs). A number of techniques have been developed for approximating number-resolving detection using these 
inherently non- number-resolving detectors [T3, EH, E2, EH EH EH EE EH EH EH j ah of which are variations on the 
idea that the optical field is distributed across multiple modes which are measured independently. The experiment is 
arranged so that the probability of any given mode being populated by more than one photon approaches zero. Thus 
in the absence of dark counts, afterpulsing or imperfect quantum efficiency, the sum of the detection events across 
the modes closely approximates the number of photons in the incident state. 

In the case of the ideal TV-port multiplexer with ideal binary photon detectors, the maximum number of photons 
able to be resolved is equal to the number of output modes. This limit arises due to excess input photons resulting in 
> 1 photons in the output modes, which are not distinguishable from the single photon case. In a realistic detector, 
however, the presence of loss means that n input photons may result in m clicks where n > m. If m < N and 
n > N then it can be seen that the loss gives rise to a useable increase in the detector's dynamic range. In general 
though, the accuracy of the reconstruction degrades and the statistical errors increase with increasing amounts of 
loss. Similarly, fragile non-classical effects are washed out by increased losses, as quantified in (2^]. This situation is 
further complicated by the presence of dark counts and afterpulsing. 

There have been many demonstrated solutions to the reconstruction of the photostatistics (photon number dis- 
tribution) of an input field using binary (on/off) detectors, but they may each be loosely grouped into one of three 
categories. All realistic photodetectors exhibit loss, which is modelled by preceding an otherwise ideal detector with a 
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beamsplitter of transmission 77. The photostatistics of the input state are thus convolved with the Bernoulli (binomial) 
distribution arising from the stochastic photon selection process performed at the beamsplitter [23|. 

This relationship forms the basis of the first category of reconstruction techniques - those involving inversion of the 
Bernoulli distribution. The analytic inversion is cumbersome and requires r\ > 0.5 [24[, an impossible constraint on 
N-ports with N > 2 which is further exaggerated at telecommunications wavelengths due to the very small quantum 
efficiencies. Numerical approaches such as maximum likelihood techniques permit more realistic values of r\ and have 
been successfully applied in [25|, [26| . 

The second category involves reconstruction of the input density matrix, p, via conditional probabilities and appears 
the most common contemporary approach. An N-port detector may be completely described by the corresponding 
conditional probability matrix P(m\n) which quantifies the probability of observing to clicks, given n incident photons. 
The experiments of [Til [l6l [27j perform a large number of measurements of an unknown input state, and record the 
number of clicks obtained from each measurement. By applying either the inverted matrix of P(m\n) values to the 
observed "number of click" statistics or utilising maximum likelihood techniques, the input photon number distribution 
may be obtained. 

Finally, as observed by [28] , an improvement in the accuracy of the reconstruction is possible by noting the coin- 
cidences of the to detection events of the N detectors. This is equivalent to recording detection signatures for each 
measurement and forms the basis of our reconstruction implementation. To date, experimental work in this third 
category has considered only N = 2 distinct output modes [291 ] for a number of detection efficiencies 77. The purpose 
of this paper is to present experimental results for N — 9 time-domain modes with a fixed (small) value of 77 for each 
mode. 



II. NUMERICAL METHODS 



Irrespective of the reconstruction category, the three approaches share the same goal common to all tomography 
problems - the calculation of the input state from a large data set of marginal distributions. Although the analytic 
relationship between the measured marginal probability distributions is usually a known function of the input state, 
this is not the case for the inverse in general. Even if an inverse analytic solution exists it is generally of no use as 
experimental and measurement uncertainties create an overcomplete data set for which there is no unique solution for 
p. The problem may thus be restated as the determination of the best estimate p for p given the measured data set. 

The estimation of p is significantly simplified by a priori assumptions of aspects of the statistics of the input 
state. If it may be assumed that the input state consists of negligible contributions from higher order Fock states, i.e. 
(k\p\k) w OVfc > K, then the reconstructed Hilbert space may be restricted to dimension (K+l) x (K+l). Additional 
constraints such as the expected distribution (e.g. Poissonian) and physicality via Tr[p] = 1 and pkk > V k < K 
further aid the reconstruction process. 

Detection with on/off photodetectors will result in a particular pattern of detector clicks, which we call a signature. 
The probability, P s i g (c?), of the detection signature occurring whereby each of the detectors in vector d trigger, and 
all others do not is given by [3(| 
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where p c (i) represents the probability of an input photon being coupled out of the i th port. Similarly, p\ oss (j) 
indicates the probability of photon loss associated with port j. Both of these parameters are a function of the 
detector architecture and are detailed in (30j . The dark count probability p^c depends upon the particular on/off 
detector employed. 

The net probability of detecting to photons, given by summing over all combinations of detection signatures where 
\d\ = to, for \d\ is total number of clicks, is thus 

P(m\n) = J2 P *iM (2) 

\d\—m 
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It is worth noting that P s i g (d) is not necessarily the same for all d where \d\ = m. That is, different signatures 
corresponding to the same measured number of photons needn't have equal probabilities of occurring. 

For the three scenarios above, the measured data set is reduced to a vector of observed probabilities for a detection 
in the relevant basis. In the binomial reconstruction approach, this corresponds to the probability of a click p c i; c k for 
each of the experimental values of detection efficiency r/x . The probability p out (m) of obtaining m output photons 
depends on the m th and subsequent terms of the input state density matrix p. i.e. 



(m) 
Pout 



(3) 



Incorporating the independent dark count statistics and assuming a binary on/off detector (typically an SPDM) 
with a linear response to photon number, Eq. [3] may be integrated over the truncated Hilbert space to give 



fc fc 



Pclick(?7, P) = Pdc + (1 - Pdc) ( „ J 7 ^ 1 - VY^Pn 

n—1 j>n 



(4) 



Corresponding expressions can be obtained for the photon counting and signature approaches from Eqs. [T]and[2] 
respectively, i.e. 



(5) 
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As stated earlier, these expressions allow the calculation of the expected statistics for a known density matrix. To 
reconstruct the unknown state, the statistical distance between the predicted and measured data sets is quantified. 
A correction is thus determined and applied to the current state prediction to create a new estimate. This process is 
repeated until either the desired degree of convergence between the datasets is achieved or the incremental changes 
to the estimated state fall below an arbitrary threshold. 

The determination of the optimal means of quantifying the statistical distance between the data sets is an area of 
ongoing research. Similarly, the calculation of the corresponding amendment to the estimated state is also a non- 
trivial exercise. A comparison of the popular contemporary methods of maximum likelihood and maximum entropy 
approaches are discussed in the context of the inversion of photostatistics data in [3l| . 

For our purposes, we quantify the distance between the theoretical and measured data sets via the mean square 
error between the distributions used as a measure of convergence by [32[ . If p and p are the calculated and observed 
marginal probability distributions in the basis applicable to the reconstruction category, and <S* is the total data set 
size, 



s 

e = J2(Pi-Pif ( 7 ) 

i 

The problem may hence be expressed as a multidimensional numerical minimisation problem, where the goal is the 
minimisation of e by appropriate selection of the (K +1) diagonal coefficients of p. For this task, the simplex method 
of Nelder and Mead |33| was chosen for ease of implementation [34[ and the requirement for evaluation of e only, and 
not its derivatives. Although convergence of the algorithm to a global minimum is not guaranteed 35], in practice its 
convergence was found to be fast and entirely adequate for a wide range of expressions for e. 
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III. 1550 NM LOOP DETECTOR 



The TDM loop detector architecture of [15j was selected because of its inherent flexibility in permitting a variable 
number N of output modes. The architecture is shown in Fig. [TJ 
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FIG. 1: Experimental configuration of the loop detector and diode laser source of pulsed coherent states. The fibre loop may 
be broken (disconnected) at point 'X' for characterisation purposes as described in the text. 



The experiment allows the various reconstruction methods to be assessed against the Poissonian photostatistics of 
the coherent output of the diode laser. The laser is biased below threshold, and pulsed briefly into operation by the 
application of a variable amplitude 20 ns pulse. 95% of the laser output is directed to a NIST-traceable Newport 3227 
power meter for calibration/renormalisation purposes. The remaining 5% is then transferred via a scries of calibrated 
fibre attenuators providing ~ 70 dB of attenuation into the loop detector. 

The loop detector comprises the input switch, output coupler and fibre delay loop. The switch was provided by 
Civcom Inc. and provides a remarkably low insertion loss (1.2 dB) combined with fast (270 ns) switching speed. As 
the fibre provides 10.044 fis of delay, the switch thus has adequate time to settle in the alternative state (shown by 
the dotted line of Fig. [1]) after admitting the pulse to the detector. The loop consists of 2 km of Corning SMF-28 
fibre, exhibiting only 0.8 dB loss at 1550 nm, with each coupler exhibiting 0.5dB of loss. The long delay allows for 
the recovery of the APD used in the SPDM following a gate pulse, in order to minimise the detrimental effects of 
afterpulsing. A portion of the beam is tapped off from the circulating field via the fibre coupler at each pass and 
directed onto an id-Quantique id-200 SPDM operating with a 20 ns gate window and a 100 kHz clock rate. This 
detector had a quantum efficiency of 10 % and a dark count probability per gate pdc = 9.6 x 10~ 4 . 

N = 14 bins was selected with a coupling coefficient p c — 0.1, yielding a possible truncation error of 23% [30] . In 
practice, this error is not observed due to the excess loop losses preventing photons making it into the last few bins 
for the range of input intensities employed. Combined with the r]d = 10% quantum efficiency of the id-Quantique 
id-200 SPDM used, the equivalent detection efficiencies Tji for each time bin are shown in Table [U 



TABLE I: Detection efficiencies corresponding to loop TDM bins. 



Bin 
V 


1 

6.83e-3 


2 

3.58e-3 


3 

1.88e-3 


4 

9.91e-4 


5 

5.22e-4 


6 

2.75e-4 


7 

1.45e-4 


Bin 
V 


8 

7.61e-5 


9 

4.01e-5 


10 
2.11e-5 


11 
l.lle-5 


12 
5.84e-6 


13 
3.08e-6 


14 
1.62e-6 



Due to the long delay employed, the loop path length is much greater than the coherence length of the laser. 
Consequently, the tapped field is uncorrelated with the circulating component within the loop, and any quantum 
correlations between successive time bins are destroyed. This allows the loop data be considered mathematically 
classical and photons to be treated as discrete particles. 

Accurate timing of the experiment is critical to ensure the APD gate window that defines each time bin is correctly 
aligned with the pulse circulating within the loop. Control of the timing is achieved via a synchronous finite state 
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machine (FSM) constructed within a field programmable gate array (FPGA). A measurement cycle consisting of the 
generation of an optical pulse, recording 14 passes of the loop and transferring the data to file for post processing 
requires 376/^s, allowing 2660 distinct measurements per second to be performed. 

a) b) 




14 15 16 17 18 19 20 21 22 23 24 25 99.52 99.54 99.56 99.58 99.6 99.62 99.64 99.66 99.68 

SPDM trigger delay (ns) Clock frequency (MHz) 



FIG. 2: Determination of a) optimum SPDM trigger delay and b) clock frequency for the loop detector. The SPDM trigger 
delay is chosen to maximise the detection probability of the first time-bin. The clock frequency is chosen to maximise X^i=2P»- 

Compensation for the multitude of electro-optic delays within the experiment is achieved by adjusting the triggering 
time of the SPDM and the clock frequency of the FSM. Fig. [2] indicates the frequency of detection events in the 
first (and subsequent) time bins as a function of time delay and clock frequency respectively. The final timing 
parameters chosen thus allow the detection efficiencies to approach those listed in Table fl] Fig. [3] illustrates the 
resulting exponentially decaying click probabilities obtained for each bin. 
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FIG. 3: Raw time bin detection (click) probabilities for r\d = 10%, n = 8.1 photons/pulse. 



In practice, it was noted that very few valid detection events were observed in the later time bins. This is as 
expected due to the vanishingly small effective detection efficiencies. Consequently, the analysis was terminated at 
N = 9 time bins. To compensate for the induced truncation error, a non-experimentally observable 10 th bin was 
re-introduced to the expression of Eq. [5] to serve as a "catch all" for photons normally lost to the analysis. 

To accurately model the 9 time-bin experimental data, Eq. [5] is thus rewritten 



Psi s (d, p) N =9 «Psi g ([d 0},p) N=10 +p sig ([d l],p) N=10 



(8) 
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and the 10 th element of the coupling and loss vectors are such that p c (10) = 1 — Yli=iPc{i) an d Pioss(lO) = 0. 
As a test of validity, this ensures 

^2p sig (d, p) N =9 = 1 V n, p (9) 

d 



IV. CHARACTERISATION 

In addition to the bin efficiencies listed in Table [H our detector can be completely characterised by its conditional 
probability matrix. The first few terms are listed in Table [TT1 and plotted graphically in Fig.UJ 



TABLE II: P(m\n) for loop detector. 



m \ n 





1 


2 


3 


4 


5 





0.9913 


0.9776 


0.9640 


0.9506 


0.9374 


0.9244 


1 


0.0086 


0.0222 


0.0356 


0.0487 


0.0614 


0.0738 


2 


3.29e-5 


1.39e-4 


3.70e-4 


7.24e-4 


1.19e-3 


1.78e-3 


3 


7.39e-8 


4.29e-7 


1.64e-6 


4.45e-6 


9.59e-6 


1.77e-5 


4 


1.06e-10 


7.90e-10 


3.95e-9 


1.39e-8 


3.79e-8 


8.56e-8 


5 


1.02e-13 


9.24e-13 


5.71e-12 


2.50e-ll 


8.30e-ll 


2.23e-10 




m 

n 



FIG. 4: Characteristic conditional probabilities P(m\n) for the 1550 nm loop detector. Note the logarithmic vertical axis used 
to accommodate the rapidly diminishing higher order terms. 

Immediately apparent is the departure from the ideal P(m\n) = S m . n , arising as a direct consequence of the loop 
losses and the poor InGaAs APD efficiency. The quantum efficiency of the id-200 can be doubled at the cost of an 
order of magnitude increase in dark counts, i.e. at rjd = 10%, pdc = 9.66 x 1CP 4 , while at rjd — 25%, pdc = 2.75 x 10~ 2 . 
However, dark counts adversely affect fidelity to a greater extent than reduced quantum efficiency |30j. Consequently, 
the lower of the id-200 detection efficiencies were used. 

The remnant dark count "noise floor" is visible at the exponential tail of Fig. [3] count probabilities. Fig. [5] provides 
further insight into the limitations of the photodetection process, illustrating the effect of disconnecting the loop at 
point X as shown in Fig. [T] for rjd — 10% and rjd — 25% respectively. The non-zero detection probabilities in bins > 2 
are the result of dark counts and after pulsing. 
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a) b) 
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FIG. 5: Measured time-bin detection probabilities for a) rjd = 10% and b) rjd = 25%, with the loop open-circuit. 

Afterpulsing describes the mechanism by which charges trapped within defects of the APD semiconductor give rise 
to spontaneous avalanche events during the next gate pulse, even in the absence of a detected photon. Unfortunately, 
since avalanches induced by a valid photon detection event are indistinguishable from APD defect/intrinsic dark count 
avalanche events, afterpulsing can artificially inflate the effective dark noise floor of many subsequent detector time 
bins. For fixed APD temperature and gating period, the afterpulsing probability scales non-linear ly with rjd [36]. This 
effect is clearly shown in Fig.(SJb), where photons are detected only during time bin 1 and afterpulsing is responsible 
for the exponential decay to the dark noise background, which is substantially greater than the rjd = 10% case of Fig. 
Eta). 



0.8 




23456789 10 
Bin number 



FIG. 6: Measurement of afterpulsing from time series data of Fig. \5\(rjd = 25%). 

In order to quantify the effects of afterpulsing, the asymptotic background value dark noise was assumed to be 
independent of detection events and was subtracted away from the open loop data. Fig. [5] indicates the result 
of comparing the probability of a click in each bin to the one before for the data of Fig. [DJb) . A significant 60% of 
detection events in subsequent time bins are thus observed to result solely from afterpulsing effects. The approximately 
constant envelope of probabilities suggests the simple "regenerative afterpulsing" model described in [37J adequately 
models our experiment. Assuming pdc <C 1, the detected p c iick data in bin i may thus be corrected to first order by 
subtracting a portion of the counts p a detected in a prior time bin, i.e. 
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Pclick(0 



Pdc + Pciick (i ) [1 - Pdc] ; i = 1 

Pciick(i) - Pdc - PaPciick(i - 1) ; i > 2 



Note that this correction aids reconstruction via the binomial method only, as integration of the comprehensive 
afterpulsing model of [37| into the expressions of Eqs. [5] and [5] is non-trivial. 



a) 



b) 
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FIG. 7: Raw and corrected time bin detection probabilities for a) r\d = 10% and b) rjd = 25% with n = 3.25 photons/pulse. 

Fig. E| shows the corrections of Eq. IIVI applied to a typical p c iick data set. Note that bin 2 consistently experiences 
the greatest correction as a consequence of afterpulsing following the significant number of detection events in the 
first time bin. With rjd = 10%, p a was measured at a modest ~ 3%. Given our inability to presently account for 
afterpulsing in the signature based reconstruction analysis, rjd = 10% was chosen for all subsequent experimental data 
runs. 



V. ANALYSIS OF RESULTS 



In order to evaluate the proposed signature-based reconstruction method, a number of coherent state input tests 
were performed. 1.5 x 10 5 fixed amplitude pulses were acquired by the loop detector for each of a range of calibrated 
input intensities. The data was then processed to form two probability distributions - p click (i) for time bins (and thus 
effective 77) i = 1 ... 10 and p S i g (d)N=9 for all 2 N = 512 signatures. 

The photostatistics of the input states were assumed to be Poissonian. For each input amplitude, the MSE given by 
Eq.[7]was calculated from the experimental data set and a range of predicted data sets generated by Eq.[4]and Eq.[6]for 
the binomial and signature schemes respectively. To create the predicted data sets, a spectrum of estimated density 
matricies were generated with Poisson-distributed photon numbers where n was swept over the interval [0, ~ 15] 
photons/pulse. Each of the evaluated input states were renomalised for a unity trace to account for the truncated 
Hilbert space, with the diagonal coefficients thus given by 



JT-T (10) 

2j<=1 i! 



By means of illustration, Fig. |H1 shows the predicted magnitudes of the first 30 p s i g probabilities for an input coherent 
state corresponding to n = 8 photons/pulse. These probability values were then compared to the measured experi- 
mental probabilities via Eq.[7]to determine the associated MSE. The indicated binary signature vectors correspond to 
the base-10 signature index, stored in bin-reversed order, i.e. time bin 1 is the right most bit 39]. For example, the 
signature vector corresponding to a detection in the first and third time-bins is d = {0, 0, 0, 0, 0, 0, 1, 0, 1}. This may 
be represented as a base-2 (binary) number, i.e. 0000001012, which has an equivalent base-10 (decimal) signature 
index 5io- Signature index 610 thus corresponds to to = 2 detection events - one in each of the second and third time 
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FIG. 8: Predicted 9-bit signature probabilities for r]d — 10%, n = 8 photons/pulse. 
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FIG. 9: Mean-square-error for coherent-state fits to raw, corrected and signature probability data for n = 6.5±0.8 photons/pulse. 



bins and so on. As expected from the P{m\n) data of Section HVl the probability of multiphoton detection events 
decreases rapidly with increasing to, hence the tiny value shown for index 23io = OOOOIOIII2, with to = 4. 

Fig. [9] illustrates the calculated values for e against mean photon- number n for the experimental binomial p c iick and 
p s ig signature probabilities for an input state of n = 6.5 ± 0.8 photons/pulse. A third trace illustrating the results 
of considering the p c iick data corrected for afterpulsing as per Eq. IIVI is also shown. Note that the latter trace 
exhibits a very much smaller minimum value of e. This suggests a much closer agreement between the modelled and 
detected data, and the hence the validity of the applied correction. The minimum value of e for a given reconstruction 
scheme indicates the value of n that results in optimum agreement between the predicted and experimentally observed 
detection probabilities. The offsets between the minima of the three curves result from the inherent sensitivity of the 
binomial approach to the absolute values of r\ used in Eq. |4j While approaches to self-calibration have been proposed 
to mitigate these effects (f| [HH , the signature reconstruction method appears relatively insensitive to variations of 
individual values of r\. For example, 5% variations of the switch loss term t s shifted the signature minimum by a 
corresponding amount, yet induced > 10% shifts in the binomial ]3 c ii C k minimum. Note also that the signature e value 
is a more sharply peaked function, potentially improving the accuracy of photon number estimation. 

The significance of the departures of the estimated values of n from the actual input state is examined in Fig. [TO] Fig. 
[TO] shows the photon number distributions corresponding to the values of n identified by the (afterpulsing-corrected) 
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a) b) 




5 10 15 5 10 15 
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FIG. 10: Photon- number distributions corresponding to the observed e (MSE) minima of Fig. [9] The best fit to the a) corrected 
binomial data occurs with n = 4.54 photons/pulse, and the b) signature data at n = 5.95 photons/pulse. The photostatistics 
(and derived error bounds) of the n — 6.5 ± 0.8 photons/pulse, input state is shown for comparison. 

binomial and signature approaches of Fig. [9] For each case, the input distribution is shown for comparison. The error 
bars were determined from the calibrated experimental uncertainty of the input amplitude. The severe consequences of 
the calibration-induced errors of the binomial scheme are starkly visible, with the estimated distribution significantly 
different from the input state. Conversely, the estimated signature statistics closely agree with the actual input 
statistics. 




FIG. 11: Estimated mean photon numbers for binomial probability and signature based reconstruction techniques. 



Fig. [TT] indicates the inferred photon numbers corresponding to the minima of the MSE e error term for the two 
approaches over a range of input amplitudes. The ideal 1:1 relationship is also shown. The signature approach 
consistently returns values within the error bounds of loop detector losses and laser intensity calibration. The linearly 
increasing error corresponds to a systematic calibration error of ±15% or 0.6dB, readily accounted for in the multiple 
insertion loss uncertainties arising from fibre-fibre interconnections. 

Finally, we investigate the performance of the signature reconstruction approach when no a priori assumptions are 
made regarding the form of the input photostatistics. In this case, we seek to find the density matrix p that minimises 
the MSE error e by searching (K + l)-dimension parameter space for an optimum set of diagonal coefficients. While 
we make use of the same coherent state experimental data used in the analysis above, the diagonal photon-number 
terms of the reconstructed density matrix are considered independent from one another and not bound to assume 
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a) b) 




Photon number Photon number 

FIG. 12: Reconstructed photon-number distributions for an input coherent state with a) fii = 4.6 ± 0.8 and b) 6.5 ± 0.8 
photons/pulse. Poissonian distributions representing the same mean photon number as the reconstructed states (n r = 4.3 and 
6.0 photons/pulse, respectively) are shown for comparison. 

Poissonian form. To this end, we employed the numerical multidimensional simplex minimisation algorithm described 
in Section [IT] to determine an optimal estimate of p that minimises e. The expression for MSE shown in Eq. [7] was 
modified to incorporate unity trace constraints 

e = ( Tr L°] - !) 2 + ^2(Psig{d,p)N=9 -Psigid, P)n=9) 2 (11) 

d 

in a manner that makes no other assumption about the physicality of the estimated state. The convergence properties 
of this expression in this application have not been comprehensively inv estig ated. Similarly, Eq. [IT] makes no attempt 
at optimisation of the likelihood or entropy of the reconstructed state [31| ■ This is identified as an area for further 
investigation. Nevertheless, the results of reconstruction of the input state photostatistics with a) fii = 4.6 ± 0.8 and 
b) 6.5 ±0.8 photons/pulse, are quite impressive, as shown in Fig. [T2l 

To permit validation of the reconstructed distributions, the mean photon numbers of the reconstructed density 
matrices were determined via n — ^2n = Qnp nn . Poissonian distributions with the corresponding mean photon num- 
bers are shown beside the reconstructed distributions of Fig. [T^for comparative purposes. The agreement between 
the coefficients is remarkably close, confirming the practical suitability of the proposed signature approach for the 
reconstruction of photostatistics. 

VI. CONCLUSION 

We have the described the use of photon count signatures from an N-port photon counting device for the purposes 
of reconstructing the photon number probability distribution of an input state, and demonstrated this technique 
experimentally on a 9-port system. We have shown that exploiting signature probabilities offers superior performance 
to approaches reconstructing photon number probabilities from the measured binomial statistics. Our method offers 
advantages in reduced statistical uncertainty, relative insensitivity to calibration errors and reduced reliance on a-priori 
assumptions about the input state. 
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